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ABSTRACT 

We present a Python package LDTk that automates the calculation of custom stellar limb 
darkening (LD) profiles and model-specific limb darkening coefficients (LDC) using the li¬ 
brary of PHOENlX-generated specific intensity spectra by |Husser et al. ( 2013) . The aim of the 
package is to facilitate analyses requiring custom generated limb darkening profiles, such as 
the studies of exoplanet transits-especially transmission spectroscopy, where the transit mod¬ 
elling is carried out for custom narrow passbands-eclipsing binaries (EBs), interferometry, 
and microlensing events. First, LDTk can be used to compute custom limb darkening pro¬ 
files with uncertainties propagated from the uncertainties in the stellar parameter estimates. 
Second, LDTk can be used to estimate the limb-darkening-model specific coefficients with 
uncertainties for the most common limb-darkening models. Third, LDTk can be directly in¬ 
tegrated into the log posterior computation of any pre-existing modelling code with minimal 
modifications. The last approach can be used to constrain the LD model parameter space di¬ 
rectly by the LD profile, allowing for the marginalization over the LD parameter space without 
the need to approximate the constraint from the LD profile using a prior. 


Key words: Methods: numerical-Planets and satellites: general-Gravitational lensing: 
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1 INTRODUCTION 

Stellar limb darkening (LD) is a major source of uncertainty in the 
characterization of transiting exoplanets (Csizmadia et ah 20 f3l|Es-] 
pinoza & Jordan 2015 1 , interferometry-based stellar radius estima¬ 
tion i White et al.|2013| Neilson & Lester|20 1 3af , and analysis of 
microlensing events (An et al.|2002) . 

Limb darkening is approximated by a variety of limb darken¬ 
ing models, or laws (Mand ef& Agol|2002||Gimenez|2006l ). These 
models aim to reproduce the stellar intensity profile as a function 
of /Qwith a relatively small number of parameters, limb darkening 
coefficients. In the context of exoplanet transit modelling, the early 
transiting exoplanet studies generally fixed the limb darkening co¬ 
efficients to values derived from the model and passband-specific 
fits to numerical stellar models by |Claret| ( |2000| |2004) . However, 
fixing the coefficients will introduce biases in the parameter esti¬ 
mates if the models fail to reproduce the true stellar intensity pro¬ 
files ( jCsizmadia et al.j|2013j [Espinoza & Jordan||2015fr , which has 
been shown to be the case in several situations where inferences 
about the true limb darkening profile have been possible JFields 


et al. |2003|[CTare t 2008 2009, ll owarth|20 lT| but see also [Miiller 


et al.|2013). In addition, even with an accurate model, fixing the 


uncertainties in the stellar parameter estimates are not propagated 
into uncertainties in the limb darkening profiles. 

The problems with fixed limb darkening have led to the use 
of more robust approaches to transit modelling, and currently the 
limb darkening is usually either constrained using (more or less 
informative) priors based on the tabulated limb darkening coef¬ 
ficients, or using uninformative priors and marginalizing over all 


the limb darkening profiles allowed by the observations (Kipping 


(2013 

2014), 

Sing (2010 

', Neilson & Lesterj(2013a b), 

Magic et al. 

(2014 

I, and f 

usser et al. 

(2013) are based on more advanced stel- 


Claret et al. 


lar atmosphere models than the early ones, using a spherical model 
geometry instead of a plane-parallel one i Claret et al.|2012]|2014| 


coefficients would lead to underestimated uncertainties since the 


|Neilson & Lester|2013a|b| |Husser et al.||20 1 3 }, usi ng fully three- 

dimensional models geometries ( |Hayek et al.|20~i~2) , and including 
hydrodynamical s imulations to account for the stellar granulation 
I Magic et al. 2014} r] 

Allowing the data to speak for itself-using noninformative pri¬ 
ors and marginalizing over the limb darkening allowed by the data- 
yields the most reliable parameter estimates, but is not the optimal 
approach in some situations. For example, in transmission spec¬ 
troscopy, where minute changes in the planetary radius as a func¬ 
tion of wavelength are used to make inferences about the planet’s 


*hannu.parviainen@physics. ox. ac.uk 
1 Where // = \/\ — z 2 = cos 7 , z is the normalized distance from the 
centre of the stellar disk, and 7 is the foreshortening angle. 


2 Note that even with reliable stellar models, errors can arise from the limb 
darkening model fitting. See Sect. 2.2 in |Espinoza & Jordan|i20 l~5j , espe¬ 
cially the need to renormalise the radius from the spherical models. 
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atmosphere, the uncertainties from unconstrained limb darkening 
easily drown any variation in the radius estimates. Further, the use 
of pretabulated coefficients is impossible due to the use of nar¬ 
row non-standard passbands, and unaccounted-for small-scale vari¬ 
ations in the limb darkening can lead to systematic trends and spu¬ 
rious features in the inferred transmission spectrum. 

A more robust way to account for limb darkening in any mod¬ 
elling situation is thus to generate limb darkening profiles from the 
modelled stellar spectra on a case-by-case basis for each passband 
in a way that propagates the uncertainties in the stellar properties 
to uncertainties in the limb darkening profile, and allows us to set a 
factor telling how much we trust on the stellar atmosphere models 
used to create the spectra. This approach offers a compromise be¬ 
tween tightly constrained and completely unconstrained limb dark¬ 
ening. The limb darkening profiles can be used to either create prior 
distributions for the limb darkening model coefficients by estimat¬ 
ing the coefficient posteriors; the profiles can be used to directly 
constrain the limb darkening in the modelling process, fitting a limb 
darkening model to the profile simultaneously with the rest of the 
analysis; or the profiles can be used directly as is, if the modelling 
approach allows for arbitrary stellar limb darkening profiles. 

We present a Python package LDTk that automates the calcu¬ 
lation of custom stellar limb darkening profiles and model-specific 
limb darkening coefficients from the library of PHOENIX-generated 
specific intensity spectra by Husser et al.|( |2013| l. LDTk propagates 
the uncertainties in the stellar parameter estimates into the uncer¬ 
tainties in the limb darkening profile, and offers methods to create 
multivariate normal priors for the most common limb darkening 
laws. The package can also be directly integrated into the mod¬ 
elling code, bypassing the need to pre-calculate the priors. While 
the following discussion considers mainly exoplanet transit mod¬ 
elling, the package can be equally used in the context of eclipsing 
binary, microlensing, and interferometry studies, or in any work 
benefiting from custom-built stellar intensity profiles. 

2 EXAMPLES 

2.1 Calculation of model coefficients 

We start with an example of limb darkening model coefficient esti¬ 
mation. LDTk offers methods to calculate the maximum likelihood 
estimates-as well as the full posterior distributions-for the coeffi¬ 
cients of the limb darkening models listed in Table|7] 

Limb darkening profiles and quadratic model coefficients for 
three simple passbands for a star with Te® = 6400 ± 100 K, 
logg = 4.5 ± 0.1, and z = 0.25 ± 0.05 can be calculated as 

from ldtk import (LDPSetCreator, 

BoxcarFilter) 

filters = [BoxcarFilter('a',450,550), 

BoxcarFilter('b',650,750), 

BoxcarFilter ('c', 850, 950)] 

sc = LDPSetCreator(filters, 

teff=[6400, 100], 

logg=[4.50, 0.10], 
z=[0.25, 0.05]) 

ps = sc.create_profiles(nsamples=500) 
qc,qe = ps.coeffs_qd(do_mc=True) 


We first import the LDPSetCreator class that will create a set 


of limb darkening profiles for a given set of filters. Next, we define 
the filters, here using simple boxcar filters that have zero trans¬ 
mission outside the given minimum and maximum wavelength 
range. We continue by creating an instance of LDPSetCreator, 
initialising it with the filter set and stellar parameter estimates. 
LDPSetCreator then creates a list of needed spectrum files, 
checks them against the local cache, and downloads the un¬ 
cached files from the Husser et al. ( |2013| > FTP library. After this, 
sc. create_prof iles is used to create the limb darkening 
profiles for each filter, shown in Fig. |T| contained in an instance 
of LDP Set class. The LDP Set class is finally used to estimate 
the quadratic limb darkening model coefficients qc and their un¬ 
certainty estimates qe, illustrated in Fig. [2] Figure [3] shows the 
quadratic model coefficients for the same star, but calculated for 
19 narrow (15 nm wide) passbands from 500 to 800 nm. 

2.2 Integration into transit modelling 

Rather than calculating the limb darkening coefficients prior to 
modelling. LDTk can be integrated directly into the code carry¬ 
ing out the log posterior evaluation. This can be done in two ways: 
(1) LDTk can be used to evaluate the log likelihood for the limb 
darkening model given model coefficients during the posterior sam¬ 
pling; (2) LDTk can be used to create a multivariate normal prior 
automatically in the initialisation before the sampling. 

2.2.1 Log likelihood evaluation 

The log likelihood for an LD profile can be calculated as 

lnl_ld = ps.lnlike_xx(cf) 

where cf is a coefficient array, and xx is the model abbreviation 
from Table[I] The log likelihood call can be integrated into the log 
posterior computation directly 

class LPFunction(object): 

def _ _init (self, ldfile=None, ...): 

« OMITTED » 

<< set up the data, models, >> 

<< parameterisation, etc. >> 

if not exists(ldfile): 
sc = LDPSetCreator(...) 
self.ps = sc.create_profiles() 
self.ps.save(ldfile) 

else: 

self.ps = load_ps(ldfile) 
self.lnlike_ld = self.ps.lnlike_qd 

def lnprior(self, pv): 

<< calculate the prior >> 

def lnlike(self, pv): 

<< calculate the log likelihood >> 

<< for the main dataset(s) >> 

def _.call._(self, pv) : 

return (self.lnprior(pv) + 
self.lnlike(pv) + 
self.lnlike_ld(pv[ldsl] ) ) 

where LPFunction is a callable class that encapsulates the log 
posterior computation. The limb darkening profiles are created dur¬ 
ing the first initialisation, and saved to a file so that we do not need 
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Figure 1. Three limb darkening profiles created by the example code in Sect. |2. 1 1 as a function of n (left) and 2 : (right). The large dots show the profile 
median with the original (Husser et al,|2013) sampling, the small dots show the same resampled to be linear in z, and the shaded areas show the 3cr profile 
uncertainties. The fluxes are normalized to unity in the centre of the stellar disk. 
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Figure 2. Bivariate normal density approximation for the quadratic model coefficient likelihood (showing 50, 95, and 99 percentile boundaries) and samples 
from the likelihood calculated using MCMC for the a passband from Sect. |24]fleft). The coefficient estimates with their uncertainties (99% central likelihood 
intervals) for all three passbands (right). The u errorbars are too small to be visible. 


to regenerate them every time. After the initialisation, all that is 
required is to add a call to calculate the log likelihood to the cal¬ 
culation of the log posterior. The example implements this in the 

_call_method, where pv is the model parameter vector, and 

ldsl is a slice selecting the limb darkening coefficients from the 
parameter vector. 

2.2.2 Multivariate normal prior generation 

The second approach for the direct integration of LDTk is to use it 
to construct a multivariate normal prior in the initialisation 

<< class initialization >> 
self.lnp_ld = self.ps.create_prior_xx() 

<< log posterior evaluation >> 


return (self.Inprior(pv) + 

self.lnp_ld(pv[Ids]) + 
self.lnlike_lc(pv)) 


where lnp_ld is now a function that returns the joint log density 
of multivariate normal priors estimated from the profiles. 

3 IMPLEMENTATION 
3.1 Overview 

The package aims to be simple to install and integrate into pre¬ 
existing Python-based modelling code, and depends only on stan¬ 
dard Python modules, NumPy, SciPy, and PyFITS. The code with 
IPython-notebook-based examples is freely available from 
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Figure 3. Quadratic limb darkening coefficients for 15 nm wide bins from 
500 to 800 nm. Note that the scales are different for the two figures. 


https://github.com/hpparvi/ldtk 


Table 1. Implemented limb darkening models, their abbreviations and defi¬ 
nitions. 


Name 

Abbr. 

Model 

Linear 

In 

i - C(1 - /*) 

Quadratic 

qd 

1 - Cu{l - fi) - c„(l - m) 2 

Nonlinear 

nl 

ic41-^ /2 ) 

General 

ge 


mean of the profile set, and the profile uncertainty by its standard 
deviation. 


3.3 Limb darkening coefficient estimation 

The LDPSet class offers methods of form coef f s_xx, where xx 
is a limb darkening model abbreviation from Table [T] to estimate 
the limb darkening model coefficients with their uncertainties. For 
example, the quadratic model coefficients can be obtained as 

qc,qe = ps.coeffs_qd() 


The repository contains also examples of the package’s basic usage, 
ways to integrate the package into a transit light curve analysis, and 
notebooks detailing the implementation of the code. 


3.2 Profile creation 

The limb darkening profiles are created using the PHOENIX- 
calculated specific intensity spectra by |Husser et al.| |2013j >. The 
spectra are stored in an FTP server as FITS files-a single file for 
a single set of stellar parameters-where each file contains a model 
spectrum spanning the wavelength range from 50 to 2600 nm with 
a 0.1 nm resolution for 78 values of /r. Each specific intensity file 
is ~15 MB, and the whole library has a size of hundreds of GBs. 
Flowever, downloading the whole library is unnecessary, and only 
a small subset of spectra are required for any given star. 

The profile creation starts with the calculation of the neces¬ 
sary spectrum files, continues with a check testing which of the 
required files are already in a local cache directory, after which 
the missing files are downloaded from the library FTP server. The 
LDPSetCreator class carries out these tasks during its initiali¬ 
sation, including all stellar spectra inside a 2n<r-wide cube on the 
three stellar parameters, where a are the stellar parameter uncer¬ 
tainties, and n is a freely set factor defaulting to three. 

Next, limb darkening profiles are calculated for each filter and 
stellar parameter set as 

1(d) — J QWF(X)I(fj,, A) dA, (1) 

where Q is the detector quantum efficiency, F is the filter transmis¬ 
sion, and I (/r, A) is the stellar spectrum. The profiles are renormal¬ 
ized in 2 (see |Espinoza & Jordan|2015] Sect. 2.2), and then used to 
construct a (Teh, log g, z) interpolation grid for each filter. 

The final limb darkening profiles are calculated using Monte 
Carlo sampling. The create_profiles method generates N 
stellar parameter samples from a multivariate normal distribution 
with the means and variances defined by the stellar parameter esti¬ 
mates and their uncertainties (we assume symmetrical normal un¬ 
certainty distributions), and calculates the limb darkening profiles 
for each sample. The final profile for each filter is defined by the 


where qc contains the coefficient estimates and qe their uncertain¬ 
ties for each filter. By default, the estimates correspond to the max¬ 
imum likelihood estimates, and the uncertainties to the 68% central 
likelihood intervals assuming multivariate normal likelihood with 
zero correlation between the coefficients. That is, we estimate the 
likelihood uncertainty as 


a = \ - 


d 2 lnP(a;o) 


d 2 * 


( 2 ) 


where P is the likelihood density, xo is the maximum likelihood 
estimate for x, a its standard deviation, and the derivative is calcu¬ 
lated numerically. 

More robust uncertainty estimates can be obtained using 
Markov Chain Monte Carlo (MCMC) sampling 


qc,qe = ps.coeffs_qd(do_mc=True) 


which also allows for the estimation of the covariance matrices 


qc,qm = ps.coeffs_qd(do_mc=True, 

return_cm=True) 


which can be useful when the estimates are used as priors in the 
modelling process. 

The methods also allow the MCMC parameter estimation to 
be fine tuned by setting the number of iterations, the chain thinning 
factor, and the number of burn-in iterations, and the MCMC chain 
is stored inside the LDPSet instance to allow for a more detailed 
convergence analysis. 


3.4 Log likelihood evaluation 

The LDPSet class offers the methods lnlike_xx, where xx is 
again the model abbreviation, for the evaluation of the model like¬ 
lihoods given a set of model coefficients. The joint likelihood for a 
quadratic model can be calculated as 

11 = ps.lnlike_qd([[0.35,0.15] , 

[0.23,0.11], 

[0.18,0.10] ] ) 




























LDTk: Limb Darkening Toolkit 5 


where we evaluate the log likelihood with two pairs of coefficients, 
one pair for each filter in the set considered in the first example. 

The likelihoods are calculated assuming that the profile values 
for a given are normally distributed, which leads to the usual log 
likelihood equation 


logL = 


N log 2-7T 
2 


N 

log tCJi 

i= 1 


\ ^ (Pi — Mi ) 2 

2-^ 2e 2 a 2 ’ 

X — 1 L 


(3) 


where N is the number of datapoints, Ui are the profile uncertain¬ 
ties, Pi are the profile values, M; are the model values, and e is a 
factor that can be used to increase the uncertainties based on how 
much we trust the stellar spectrum models. 


3.5 Uncertainty multiplier 

The uncertainty multiplier e is a subjective factor that defines how 
strongly the limb darkening profile (or the prior created from it) 
constrains the final analysis (that is, how much we trust the stel¬ 
lar atmosphere models used to create the profiles.) The uncertainty 
multiplier is applied to log likelihood calculations and model coef¬ 
ficient uncertainty estimation as described in Eq.[3] 


3.6 Resampling 

The sampling in g used by Husser et al.] ( |2013[ l focuses 
on the detailed sampling of the stellar edge, while the in¬ 
ner parts of the stellar disk are sampled sparsely. While sam¬ 
pling like this is a good choice for representing the stellar in¬ 
tensity profile, it may not be optimal for likelihood evalua¬ 
tion and model fitting (at least without sample weighting), since 
the fit will be dominated by the very edge. A linear sampling 
in 2, for example, may be a better approach. LDP Set includes 
a method to resample the profiles to any given vector of /.t 
or 2 values, resample (mu=mu_array, z=None), and also 
two utility methods resample_linear_mu (nsamples) and 
resample_linear_z (nsamples) to resample the profiles 
linearly in /r or z. 


3.7 Filters 

A filter defines a single passband, and can be either a function or a 
callable class that returns an array of transmission values between 
zero and one given an array of wavelengths in nm. The package 
comes with some utility filters, such as the BoxcarFilter used 
in the first example, and the TabulatedFilter, which allows 
for the use of tabulated transmission values listed as a function of 
wavelength. 


3.8 Quantum efficiency 

Any function that evaluates the quantum efficiency as a func¬ 
tion of wavelength (similar to a filter) can be used to de¬ 
fine the instrumental quantum efficiency. LDPSetCreator can 
be given a detector quantum efficiency curve in the initialisa¬ 
tion, and a TabulatedQE class can be used the same way as 
TabulatedFilter. 


4 DISCUSSION 

While our ability to model stellar limb darkening is improving due 
to the advances in the stellar atmosphere models, we still need to 
develop the ways we integrate the information from these mod¬ 
els into our analyses. Current stellar spectrum libraries allow us 
to calculate limb darkening profiles for custom passbands while 
propagating the uncertainties in the stellar parameter estimates into 
the limb darkening profile uncertainties, and we should take ad¬ 
vantage of this capability. The marginalization over all the limb 
darkening profiles allowed by the observations yields the most ro¬ 
bust (or, at least, conservative) parameter estimates, but the use of 
information from the stellar models is also justified in many sit¬ 
uations. Transmission spectroscopy, future planet-characterization 
space missions such as TESS (Ricker et al. 2014) and Plato (Rauer 
|et al.|[20T4] >, stellar interferometry, and microlensing studies will 
all benefit from more advanced ways to integrate the results from 
stellar atmosphere modelling into the analysis process. 

The LDTK-calculated profiles and model coefficients cannot 
be directly compared with the previous tabulations by Claret et al. 
< ]201 or others due to the differences in handling of the geometry 
and sampling. Especially the redefinition of the stellar edge (see |Es-| 
Ipinoza & Jordan|2015| Sect. 2.2) and the use of uniform sampling 
in 2 -space will affect the model coefficient estimates (however, 
these both can be changed by the user). Figure [4] shows LDTk- 
created profiles for a T E r e = 6500 ± 50 K, log g = 4.5 ± 0.1, 
and z = 0.0 ± 0.05 star observed in g\ r\ i', and z' filters with 
the corresponding quadratic models using the tabulations by |Claret| 
|et al.| ( |20T3) . The |Claret et al.| fits correspond to spherical PHOENIX 
models where the intensities corresponding to /r < 0.1 have been 
removed from the fit. The |Claret et al.||2013| models agree with 
the LDTk results close to the stellar limb, but deviate significantly 
for a large span of g (and thus for most part of the stellar disk) in 
the blue passbands. 


5 CONCLUSIONS 

We have presented a Python package LDTk that aims to facilitate 
the inclusion of information from the [Husser et al.[ ( |2013| l stellar 
atmosphere model spectrum library into astrophysical modelling 
problems. LDTk offers a modelling approach where the informa¬ 
tion from the stellar atmosphere models can be used to constrain 
the stellar limb darkening softly, the strength of the constraint being 
defined by the researcher’s (subjective) trust to the stellar models. 
The package can be used to calculate limb darkening coefficients 
and their uncertainties for freely defined passbands prior to the 
modelling, or it can be directly integrated into the modelling code. 
LDTk can also be used to generate limb darkening profile samples 
directly, bypassing the need to use simple limb darkening models if 
the modelling approach allows for the use of arbitrary limb darken¬ 
ing profiles (such as the BATMAN transit modelling code, |Kreidberg| 
|2015| . This can be useful when the features that are not well cap¬ 
tured by the simple limb darkening models, such as the very edge 
of the star, are important. 
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Figure 4. Limb darkening profiles for a r^fr = 6500 it 50 K, log g = 4.5 =t 0.1, and z = 0.0 ± 0.05 star observed in g', r', i\ and z 'filters (coloured areas), 
and the corresponding quadratic model fits from|Claret & Bloemen|i|20l I ji using either the least-squares (left) or flux conservation method (right). 
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